/*
 * EBFV1_loadBoundary.cpp
 *
 *  Created on: 16/09/2012
 *      Author: rogsoares
 */

#include "EBFV1_preprocessor.h"

void loadDij(ifstream &fid, TMesh *pTMesh_data){
	PetscPrintf(PETSC_COMM_WORLD,"Loading boudanry elements...\t");

	int j;
	int dim = pTMesh_data->getMeshDim();
	int pos = 2*dim + 2;

	std::string str[10];
	int ndom = pTMesh_data->getNumDomains();
	BdryElementData* edata;
	std::vector<BED_list> bdryElemVectorList(ndom);					// list for domain dom vertices

	for(int i=0; i<ndom; i++){
		BED_list list;
		bdryElemVectorList[i] = list;
	}

	getline(fid,str[0],'\n');
	getline(fid,str[0],'\n');
	while ( 1 ){
		edata = new BdryElementData;
		// flag dom1 dom2 id0 id1 id2      Dijx                Dijy             Dijz
		// 999  3300 4400 161 160 866 0.000000000000000 -0.000173423583642 -0.000000000000000
		//  0     1   2    3   4   5         6                  7                 8
		// 1000 3300 0 18 5746 0.000000000000000 -0.006578947350000
		fid>>str[0];
		if (str[0] == "end") break;
		edata->flag = atoi(str[0].c_str());
		for (j=1; j<=pos; j++){
			fid >> str[j];
		}

		edata->dom1 = atoi(str[1].c_str());
		edata->dom2 = atoi(str[2].c_str());
		edata->ID_0 = atoi(str[3].c_str());
		edata->ID_1 = atoi(str[4].c_str());
		if (dim==2){
			if (edata->ID_0 > edata->ID_1){
				std::swap(edata->ID_0,edata->ID_1);
			}
		}
		else{
			edata->ID_2 = atoi(str[5].c_str());
		}

		/*
		 * For 2-D domains, multiply Dij by reservoir height (H) for 2D/3D simulations (physics occurs only on 2-D but reservoir volume is considered)
		 */
		edata->Dij = new double[dim];
		edata->versor = new double[dim];
		double H = 1.0;//(pGCData->getMeshDim()==2)?pGCData->getReservoirHeight():1.0;
		for (j=0; j<dim; j++){
			edata->Dij[j] = H*strtod(str[j+dim+3].c_str(),0);
		}

		// before insert edata on list check to which domains edata belongs. edata can be a external boundary element or a boundary element between two domains
		// if edata is a external boundary element dom1 or dom2 must be zero, if both are different from zero, then it's located between two domains.
		int pos;
		if (edata->dom1){
			pos = pTMesh_data->getMappedIndexFromFlag(edata->dom1);
			BED_list &list = bdryElemVectorList[pos];
			list.push_back(edata);
		}
		if (edata->dom2){
			pos = pTMesh_data->getMappedIndexFromFlag(edata->dom2);
			BED_list &list = bdryElemVectorList[pos];
			list.push_back(edata);
		}
	}
	transferBdryElementDataFromListToMatrix(pTMesh_data,bdryElemVectorList);
	PetscPrintf(PETSC_COMM_WORLD,"done\n");
}

// todo: check Dij orientation: from dom1 to dom2!!!!!
void transferBdryElementDataFromListToMatrix(TMesh *pTMesh_data, std::vector<BED_list> &bdryElemVectorList){
	int ndom = pTMesh_data->getNumDomains();
	for(int i=0; i<ndom; i++){
		BED_list list = bdryElemVectorList[i];
		int nelem = (int)list.size();
		pTMesh_data->bdry_elem_list[i].allocateMemory(nelem,1);
		BED_list::iterator BED_iter = list.begin();
		int row = 0;
		for(;BED_iter!=list.end();BED_iter++){
			pTMesh_data->bdry_elem_list[i].setValue(row++,0,*BED_iter);
		}
		list.clear();
	}
	bdryElemVectorList.clear();
}

void printBdryElemData(TMesh *pTMesh_data){
	int nrows, ncols;
	int ndom = pTMesh_data->getNumDomains();
	for (int i=0; i<ndom; i++){
		pTMesh_data->bdry_elem_list[i].getSize(nrows,ncols);
		for (int j=0; j<nrows; j++){
			BdryElementData* edata = pTMesh_data->bdry_elem_list[i].getValue(j,0);
			cout << setprecision(5) << scientific;
			cout << pTMesh_data->getdomFlags(i) << "\t";
			cout << edata->ID_0 << "\t";
			cout << edata->ID_1 << "\t";
			cout << edata->ID_2 << "\t";
			cout << edata->flag  << "\t";
			cout << edata->dom1  << "\t";
			cout << edata->dom2  << "\t";
			cout << edata->numRC << "\t";
			cout << edata->Dij[0] << "\t";
			cout << edata->Dij[1] << "\t";
			cout << edata->Dij[2] << "\t";
			cout << endl;
		}
	}
	cout << "ndom = " << ndom << endl;
}
